\chapter{1767-1772年，拉格朗日点计算}


\section{拉格朗日点推导}
https://zhuanlan.zhihu.com/p/76332413
任何两个大质量天体系统，因万有引力作用而相互绕行，在它们的轨道平面上，存在着一些特殊的点。当在这些点引入小质量天体时，小质量天体受到两个大质量天体引力的合力，恰好等于它与两个大质量天体一起转动需要的向心力。这意味着，小质量天体在这些点上，能保持与两个大质量天体的相对位置始终不变。

1767年，欧拉求出了这样的3个点，记为L1，L2，L3。它们都在两大质量天体的连线上。1772年，拉格朗日又找到了另外的2个点，记为L4，L5，它们在两大质量天体连线之外，这2点关于两大质量天体的连线是对称的，它们中的任一点与两大质量天体正好构成一个等边三角形(如图\ref{LagrangePointsPos})。后来，人们习惯上把这5个点统称为“拉格朗日点”，其中拉格朗日找到的那两点又叫“三角拉格朗日点”。

\includegraphics[scale=0.5]{RocheLobeOverflow}

必须明确，位于拉格朗日点的小质量天体质量远远小于两个大质量天体，所以不考虑小质量天体的反作用对两个大质量天体运动的影响。同时，除两个大质量天体外，也不考虑其他天体对这个小质量天体的引力作用，即位于拉格朗日点的小质量天体的运动完全由两个大质量天体的引力决定。

上述求解拉格朗日点的问题，是一个典型的“平面圆型限制性三体问题”。众所周知，两个天体在相互引力作用下如何运动的所谓“二体问题”，早已获得精确而完美的解答。而“三体问题”，则是近300年来数学物理中一直未能彻底解决的难题。19世纪末，法国数学家庞加莱采用美国数学家希尔提出的简化模型：“假定有两个天体，它们在万有引力作用下，围绕共同质心，沿着椭圆形的轨道，作严格的周期性运动；另有一颗宇宙尘埃，在这两个天体的引力场中游荡。两天体可完全不必理会这颗粒产生的引力对它们之间运动的影响，因为颗粒的质量相对它们自己来说实在太小了。可是颗粒的运动会是怎样的呢?这简化模型称为‘限制性三体问题’。”即使是“限制性三体问题”，也异常复杂，庞加莱借助他天才的创造-相图、拓扑方法、微分方程定性理论，发现颗粒的运动是没完没了的自我缠结，而且高度不稳定。

“平面圆型限制性三体问题”是讨论得较多的比较简单的类型，即将“限制性三体问题”中两个大质量天体围绕质心O的椭圆轨道运动简化为圆运动。所以，拉格朗日点上的小质量天体，在两个大质量天体引力的作用下，绕质心O以同样的角速度作匀速圆周运动，这样才能与两个大质量天体保持相对位置不变。

文献[2]对“三角拉格朗日点”进行了理论验算，证明了L4、L5中任一个与两个大质量天体正好构成一个等边三角形。
\subsection{一种推导方法}
https://wenku.baidu.com/view/ed3b42ba185f312b3169a45177232f60ddcce7c5.html?rec\_flag=default\&sxts=1583254173729

考虑两个大质量$m_1\gg m_2$的星体只受到互相之间的引力作用而不受到任何其他的作用，围绕质心O作圆周运动。建立极坐标系，以质心O为极点，水平向右射线OX为极轴，逆时针方向为极角正方向。m1,m2极坐标分别是(r1,$\pi$),(r2,0)，r1,r2分别是m1、m2到质心O的距离。另外有飞行器质量$m\ll m_2$位于P(r,$\theta$)，对m1,m2引力场没有影响，受两个星体的引力作用。飞行器的拉格朗日量L为：

\begin{align}
	L&=\frac{1}{2}m(\dot{r}^2+r^2\dot{\theta}^2)+\frac{Gmm_1}{\sqrt{r^2+r_1^2+2rr_1cos(\theta-\theta_1)}}+\frac{Gmm_2}{\sqrt{r^2+r_2^2-2rr_2cos(\theta-\theta_1)}}\label{LagrangePointsEQ101}
\end{align}

继续计算：

\begin{align}
	\frac{\partialL}{\partialr}&=mr\dot{\theta}^2-\frac{Gm_1m(r+r_1cos(\theta-\theta_1))}{(r^2+r_1^2+2rr_1cos(\theta-\theta_1))^{\frac{3}{2}}}-\frac{Gm_2m(r-r_2cos(\theta-\theta_1))}{(r^2+r_2^2-2rr_2cos(\theta-\theta_1))^{\frac{3}{2}}}\label{LagrangePointsEQ102}\\
	\frac{d}{dt}\frac{\partialL}{\partial\dot{r}}&=m\ddot{r}\label{LagrangePointsEQ103}\\
	\frac{\partialL}{\partial\theta}&=\frac{Gm_1m(rr_1sin(\theta-\theta_1))}{(r^2+r_1^2+2rr_1cos(\theta-\theta_1))^{\frac{3}{2}}}-\frac{Gm_2m(rr_2sin(\theta-\theta_1))}{(r^2+r_2^2-2rr_2cos(\theta-\theta_1))^{\frac{3}{2}}}\label{LagrangePointsEQ104}\\
	\frac{d}{dt}\frac{\partialL}{\partial\dot{\theta}}&=mr^2\ddot{\theta}+2mr\dot{r}\dot{\theta}\label{LagrangePointsEQ105}
\end{align}

根据拉格朗日点的定义我们知道飞行器在拉格朗日点的时候没有径向速度，所以r对时间的求导都是0，其次两个星体作匀速圆周运动，而飞行器保持与它们的相对位置不变所以其角加速度也为0，角速度为常量等于两个星体的角速度。所以有以下方程组：

\begin{align}
	0&=mr\dot{\theta}^2-\frac{Gm_1m(r+r_1cos(\theta-\theta_1))}{(r^2+r_1^2+2rr_1cos(\theta-\theta_1))^{\frac{3}{2}}}-\frac{Gm_2m(r-r_2cos(\theta-\theta_1))}{(r^2+r_2^2-2rr_2cos(\theta-\theta_1))^{\frac{3}{2}}}\label{LagrangePointsEQ106}\\
	\frac{d}{dt}\frac{\partialL}{\partial\dot{r}}&=m\ddot{r}\label{LagrangePointsEQ107}\\
	0&=\frac{Gm_1m(rr_1sin(\theta-\theta_1))}{(r^2+r_1^2+2rr_1cos(\theta-\theta_1))^{\frac{3}{2}}}-\frac{Gm_2m(rr_2sin(\theta-\theta_1))}{(r^2+r_2^2-2rr_2cos(\theta-\theta_1))^{\frac{3}{2}}}\label{LagrangePointsEQ108}\\
	\frac{d}{dt}\frac{\partialL}{\partial\dot{\theta}}&=mr^2\ddot{\theta}+2mr\dot{r}\dot{\theta}\label{LagrangePointsEQ109}
\end{align}

两个星体的动量守恒，得到

\begin{align}
	m_1r_1&=m_2r_2\label{LagrangePointsEQ110}
\end{align}

式\ref{LagrangePointsEQ108}可以化简为

\begin{align}
	\frac{sin(\theta-\theta_1)}{(r^2+r_1^2+2rr_1cos(\theta-\theta_1))^{\frac{3}{2}}}&=\frac{sin(\theta-\theta_1)}{(r^2+r_2^2-2rr_2cos(\theta-\theta_1))^{\frac{3}{2}}}\label{LagrangePointsEQ110}
\end{align}

于是式\ref{LagrangePointsEQ110}成立就有两种情况，即分子为0且分母不为0，或者分子不为0且分母相等。
\subsubsection{计算拉格朗日点L4,L5}
式\ref{LagrangePointsEQ110}分母相等则飞行器与两个星体的距离相等，至少形成一个等腰三角形，进一步列出具体等式有：

\begin{align}
	r^2+r_1^2+2rr_1cos(\theta-\theta_1)&=r^2+r_2^2-2rr_2cos(\theta-\theta_1)\label{LagrangePointsEQ111}\\
	(r_2-r_1)(r_2+r_1)&=2r(r_2+r_1)cos(\theta-\theta_1)\\
	r_2-r_1=2rcos(\theta-\theta_1)\label{LagrangePointsEQ112}
\end{align}

由上面的定义我们知道$\theta-\theta_1$就是飞行器与星体m1的夹角，将这个r1,r2的关系代入到方程组的第一个式子中得到

\begin{align}
	Gm_1+Gm_2&=(r^2+r_1^2+2rr_1cos(\theta-\theta_1))^{\frac{3}{2}}\dot{\theta}^2\label{LagrangePointsEQ113}
\end{align}

我们又知道，m1,m2都在做圆周运动，其运动方程为：

\begin{align}
	\frac{Gm_1}{(r_1+r_2)^2}&=\dot{\theta}^2r_2\label{LagrangePointsEQ114}\\
	\frac{Gm_2}{(r_1+r_2)^2}&=\dot{\theta}^2r_1\label{LagrangePointsEQ115}
\end{align}

式\ref{LagrangePointsEQ114}加式\ref{LagrangePointsEQ115}，得到：

\begin{align}
	Gm_1+Gm_2&=\dot{\theta}^2(r_1+r_2)^3\label{LagrangePointsEQ116}
\end{align}

最终我们得到

\begin{align}
	r_1+r_2&=(r^2+r_1^2+2rr_1cos(\theta-\theta_1))^{\frac{1}{2}}\label{LagrangePointsEQ117}
\end{align}

即三个物体的连心线不仅构成等腰三角形而且是等边三角形。这样的等边三角形有两个（一上一下），所以这是两个拉格朗日点L4,L5。剩下的情况就是分子为0的情况：

\begin{align}
	sin(\theta-\theta_1)&=0\label{LagrangePointsEQ118}\\
	\theta-\theta_1&=n\pi,n=0,1,\ldots\label{LagrangePointsEQ119}
\end{align}

此时飞行器与两星体处于同一直线上，所以这种情况简单得多，可以方便地使用牛顿力学来确定具体的位置。
\subsubsection{计算拉格朗日点L1}
https://zhuanlan.zhihu.com/p/76332413
\label{CalculationL1}
计算原则：在拉格朗日点，质量受到的引力的合力=向心力。

设双星距离R=r1+r2，主星距离旋转中心r1=aR，

\begin{align}
	R&=r_1+r_2\label{LagrangePointsEQ301}\\
	r_1&=aR\label{LagrangePointsEQ302}\\
	r_2&=(1-a)R\label{LagrangePointsEQ303}\\
	a&=\frac{m_2}{m_1+m_2}\label{LagrangePointsEQ304}\\
	\frac{Gm_1}{(R-r)^2}-\frac{Gm_2}{r^2}&=\dot{\theta}^2[(R-r-aR)]\label{LagrangePointsEQ305}\\
	\frac{Gm_1}{(R-r)^2}-\frac{Gm_2}{r^2}&=\dot{\theta}^2R-\dot{\theta}^2r-\dot{\theta}^2aR\label{LagrangePointsEQ306}
\end{align}

从\ref{LagrangePointsEQ116}、\ref{LagrangePointsEQ301}、\ref{LagrangePointsEQ304}得

\begin{align}
	\dot{\theta}^2&=G\frac{m_1+m_2}{R^3}\label{LagrangePointsEQ307}\\
	\dot{\theta}^2&=G\frac{m_2}{aR^3}\label{LagrangePointsEQ308}
\end{align}

代入到\ref{LagrangePointsEQ306}得到

\begin{align}
	\frac{Gm_1}{(R-r)^2}-\frac{Gm_2}{r^2}&=G\frac{m_1+m_2}{R^3}R-G\frac{m_1+m_2}{R^3}r-G\frac{m_2}{aR^3}aR\label{LagrangePointsEQ309}\\
	\frac{Gm_1}{(R-r)^2}-\frac{Gm_2}{r^2}&=G\frac{m_1}{R^3}R-G\frac{m_1+m_2}{R^3}r\label{LagrangePointsEQ310}
\end{align}

因为$r\ll R$，所以r可以作为小量近似

\begin{align}
	\frac{1}{(R-r)^2}&\approx\frac{1}{R^2}+\frac{2r}{R^3}\label{LagrangePointsEQ311}
\end{align}

代入到\ref{LagrangePointsEQ310}得到

\begin{align}
	\frac{Gm_1}{R^2}+\frac{2Gm_1r}{R^3}-\frac{Gm_2}{r^2}&=G\frac{m_1}{R^3}R-G\frac{m_1+m_2}{R^3}r\label{LagrangePointsEQ312}\\
	\frac{2Gm_1}{R^3}-\frac{Gm_2}{r^3}&=-G\frac{m_1+m_2}{R^3}\label{LagrangePointsEQ313}\\
	\frac{3Gm_1+Gm_2}{R^3}&=\frac{Gm_2}{r^3}\label{LagrangePointsEQ314}\\
	\frac{3m_1+m_2}{R^3}&=\frac{m_2}{r^3}\label{LagrangePointsEQ315}\\	
	r&=R\sqrt[3]{\frac{m_2}{3m_1+m_2}}\label{LagrangePointsEQ316}
\end{align}

根据\ref{LagrangePointsEQ315}，略去小量，也可以改写成密度比。设星体为球形且密度分别为$\rho_1,\rho_2$

\begin{align}
	m_1&=\frac{4\pi}{3}R_1^3\rho_1\label{LagrangePointsEQ317}\\	
	m_2&=\frac{4\pi}{3}R_2^3\rho_2\label{LagrangePointsEQ318}\\
	\frac{m_1}{m_2}&=\frac{R_1^3\rho_1}{R_2^3\rho_2}\label{LagrangePointsEQ319}
\end{align}

代入到\ref{LagrangePointsEQ315}得到

\begin{align}
	3m_1/m_2+1&=\frac{R^3}{r^3}\label{LagrangePointsEQ320}\\	
	3\frac{R_1^3\rho_1}{R_2^3\rho_2}+1&=\frac{R^3}{r^3}\label{LagrangePointsEQ321}\\
	3\frac{R_1^3}{R_2^3}\frac{\rho_1}{\rho_2}+1&=\frac{R^3}{r^3}\label{LagrangePointsEQ322}\\
	3(\frac{R_1}{R_2})^3\frac{\rho_1}{\rho_2}+1&=(\frac{R}{r})^3\label{LagrangePointsEQ323}\\
	3(\frac{R_1}{R_2})^3\frac{\rho_1}{\rho_2}&\approx(\frac{R}{r})^3\label{LagrangePointsEQ324}\\
	\sqrt[3]{3\frac{\rho_1}{\rho_2}}\frac{R_1}{R_2}&\approx\frac{R}{r}\label{LagrangePointsEQ325}\\
	1.442\sqrt[3]{\frac{\rho_1}{\rho_2}}\frac{R_1}{R_2}&\approx\frac{R}{r}\label{LagrangePointsEQ326}
\end{align}

式\ref{LagrangePointsEQ320}$\sim$\ref{LagrangePointsEQ326}表示星体质量、半径、密度、轨道半径在拉格朗日点L1时需要满足的约束关系。
\subsubsection{计算拉格朗日点L2}
\label{CalculationL2}
按照\ref{CalculationL1}类似的方法，可以得到拉格朗日点L2。

因为L2点位于星体m1到m2的延长线上，在式\ref{LagrangePointsEQ305}以R+r替换R-r，并且要改变m2引力方向，于是

\begin{align}
	\frac{Gm_1}{(R+r)^2}+\frac{Gm_2}{r^2}&=\dot{\theta}^2[(R+r-aR)]\label{LagrangePointsEQ405}\\
	\frac{Gm_1}{(R+r)^2}+\frac{Gm_2}{r^2}&=\dot{\theta}^2R+\dot{\theta}^2r-\dot{\theta}^2aR\label{LagrangePointsEQ406}
\end{align}

\ref{LagrangePointsEQ307}、\ref{LagrangePointsEQ308}代入到\ref{LagrangePointsEQ406}得到

\begin{align}
	\frac{Gm_1}{(R+r)^2}+\frac{Gm_2}{r^2}&=G\frac{m_1+m_2}{R^3}R+G\frac{m_1+m_2}{R^3}r-G\frac{m_2}{aR^3}aR\label{LagrangePointsEQ409}\\
	\frac{Gm_1}{(R+r)^2}+\frac{Gm_2}{r^2}&=G\frac{m_1}{R^3}R-G\frac{m_1+m_2}{R^3}r\label{LagrangePointsEQ410}
\end{align}

因为$r\ll R$，所以r可以作为小量近似

\begin{align}
	\frac{1}{(R+r)^2}&\approx\frac{1}{R^2}-\frac{2r}{R^3}\label{LagrangePointsEQ411}
\end{align}

代入到\ref{LagrangePointsEQ410}得到

\begin{align}
	\frac{Gm_1}{R^2}-\frac{2Gm_1r}{R^3}+\frac{Gm_2}{r^2}&=G\frac{m_1}{R^3}R-G\frac{m_1+m_2}{R^3}r\label{LagrangePointsEQ412}\\
	\frac{2Gm_1}{R^3}-\frac{Gm_2}{r^3}&=G\frac{m_1+m_2}{R^3}\label{LagrangePointsEQ413}\\
	\frac{Gm_1-Gm_2}{R^3}&=\frac{Gm_2}{r^3}\label{LagrangePointsEQ414}\\
	\frac{m_1-m_2}{R^3}&=\frac{m_2}{r^3}\label{LagrangePointsEQ415}\\	
	r&=R\sqrt[3]{\frac{m_2}{m_1-m_2}}\label{LagrangePointsEQ416}
\end{align}

类似地得到

\begin{align}
	m_1/m_2-1&=\frac{R^3}{r^3}\label{LagrangePointsEQ420}\\	
	\frac{R_1^3\rho_1}{R_2^3\rho_2}-1&=\frac{R^3}{r^3}\label{LagrangePointsEQ421}\\
	\frac{R_1^3}{R_2^3}\frac{\rho_1}{\rho_2}-1&=\frac{R^3}{r^3}\label{LagrangePointsEQ422}\\
	(\frac{R_1}{R_2})^3\frac{\rho_1}{\rho_2}-1&=(\frac{R}{r})^3\label{LagrangePointsEQ423}\\
	(\frac{R_1}{R_2})^3\frac{\rho_1}{\rho_2}&\approx(\frac{R}{r})^3\label{LagrangePointsEQ424}\\
	\sqrt[3]{\frac{\rho_1}{\rho_2}}\frac{R_1}{R_2}&\approx\frac{R}{r}\label{LagrangePointsEQ426}
\end{align}

式\ref{LagrangePointsEQ420}$\sim$\ref{LagrangePointsEQ426}表示星体质量、半径、密度、轨道半径在拉格朗日点L2时需要满足的约束关系。
\subsubsection{计算拉格朗日点L3}
\label{CalculationL3}
按照\ref{CalculationL2}类似的方法，可以得到拉格朗日点L3。

这里需要思考，如果还和求L1,L2时设值一样,就不会有$r\ll R$了,难以进行小量近似求解,会很棘手。

所以我们设L3点离太阳中心$R+\delta$,这里的$\delta$便显然是个小量,因为L3点位于星体m2到m1的延长线上，距离主星m1更近，质量m3受到星体引力的合力等于向心力，于是得到

\begin{align}
	\frac{Gm_1}{(R+\delta)^2}+\frac{Gm_2}{(2R+\delta)^2}&=\dot{\theta}^2[(R+\delta+aR)]\label{LagrangePointsEQ505}\\
	\frac{Gm_1}{(R+\delta)^2}+\frac{Gm_2}{(2R+\delta)^2}&=\dot{\theta}^2R+\dot{\theta}^2\delta+\dot{\theta}^2aR\label{LagrangePointsEQ506}
\end{align}

\ref{LagrangePointsEQ307}、\ref{LagrangePointsEQ308}代入到\ref{LagrangePointsEQ506}得到

\begin{align}
	\frac{Gm_1}{(R+\delta)^2}+\frac{Gm_2}{(2R+\delta)^2}&=G\frac{m_1+m_2}{R^3}R+G\frac{m_1+m_2}{R^3}\delta+G\frac{m_2}{aR^3}aR\label{LagrangePointsEQ509}\\
	\frac{Gm_1}{(R+\delta)^2}+\frac{Gm_2}{(2R+\delta)^2}&=G\frac{m_1+2m_2}{R^2}+G\frac{m_1+m_2}{R^3}\delta\label{LagrangePointsEQ510}
\end{align}

因为$\delta\ll R$，所以$\delta$可以作为小量近似

\begin{align}
	\frac{1}{(R+\delta)^2}&\approx\frac{1}{R^2}-\frac{2\delta}{R^3}\label{LagrangePointsEQ511}\\
	\frac{1}{(2R+\delta)^2}&\approx\frac{1}{4R^2}-\frac{2\delta}{8R^3}\label{LagrangePointsEQ5110}	
\end{align}

代入到\ref{LagrangePointsEQ510}得到

\begin{align}
	\frac{Gm_1}{R^2}-\frac{2Gm_1\delta}{R^3}+\frac{Gm_2}{4R^2}-\frac{Gm_2\delta}{4R^3}&=G\frac{m_1+2m_2}{R^2}+G\frac{m_1+m_2}{R^3}\delta\label{LagrangePointsEQ512}\\
	\frac{-2m_1\delta}{R^3}+\frac{m_2}{4R^2}-\frac{m_2\delta}{4R^3}&=\frac{2m_2}{R^2}+\frac{m_1+m_2}{R^3}\delta\label{LagrangePointsEQ513}\\
	\frac{-3m_1-m_2}{R^3}\delta+\frac{-7m_2}{4R^2}-\frac{m_2\delta}{4R^3}&=0\label{LagrangePointsEQ514}\\
	\frac{-12m_1-5m_2}{4R^3}\delta-\frac{7m_2}{4R^2}&=0\label{LagrangePointsEQ515}\\				
	\delta&=-R\frac{7m_2}{12m_1+5m_2}\label{LagrangePointsEQ516}\\
	\delta&=-R\frac{7}{12m_1/m_2+5}\label{LagrangePointsEQ517}\\
	\delta&\approx-R\frac{7m_2}{12m_1}\label{LagrangePointsEQ518}\\
	R+\delta&\approxR(1-\frac{7m_2}{12m_1})\label{LagrangePointsEQ519}\\
	R+\delta&=R(1-\frac{7m_2}{12m_1+5m_2})\\
	&=R(\frac{12m_1+5m_2-7m_2}{12m_1+5m_2})\\
	R+\delta&=R(\frac{12m_1-2m_2}{12m_1+5m_2})\label{LagrangePointsEQ520}\\
	R+\delta+aR&\approxR(1-\frac{7m_2}{12m_1}+a)\\
	R+\delta+aR&=R(1-\frac{7m_2}{12m_1}+\frac{m_2}{m_1+m_2})\\
	R+\delta+aR&=R(1-\frac{7m_2}{12m_1}+\frac{12m_2}{12m_1+12m_2})\label{LagrangePointsEQ522}\\
	R+\delta+aR&\approxR(1+\frac{5m_2}{12m_1})\label{LagrangePointsEQ523}
\end{align}

也可以强行近似成公式\ref{LagrangePointsEQ518}，那么L3到主星球心的距离就是式\ref{LagrangePointsEQ519}，L3到主星和行星的旋转中心的距离就是式\ref{LagrangePointsEQ523}。

类似地得到

\begin{align}
	\frac{12m_1-2m_2}{12m_1+5m_2}&=\frac{R+\delta}{R}\label{LagrangePointsEQ520}\\	
	\frac{12R_1^3\rho_1}{7R_2^3\rho_2}&=-\frac{R}{\delta}\label{LagrangePointsEQ521}\\
	\frac{12R_1^3}{7R_2^3}\frac{\rho_1}{\rho_2}&=-\frac{R}{\delta}\label{LagrangePointsEQ522}\\
	(\frac{12R_1}{7R_2})^3\frac{\rho_1}{\rho_2}&=-\frac{R}{\delta}\label{LagrangePointsEQ523}\\
	\sqrt[3]{\frac{\rho_1}{\rho_2}\frac{\delta}{R}}\frac{12R_1}{7R_2}&=-1\label{LagrangePointsEQ526}
\end{align}

式\ref{LagrangePointsEQ520}$\sim$\ref{LagrangePointsEQ526}表示星体质量、半径、密度、轨道半径在拉格朗日点L3时需要满足的约束关系。

综上，我们遍历所有情况，得到五个拉格朗日点L1-L3,L4,L5。
\subsection{与简谐振动对比}
设飞行器与m2之间
\subsection{另外一种推导方法}
https://zhuanlan.zhihu.com/p/61552995

[Dr.S]拉格朗日点(平动点)的近似推导

Dr.Space
\subsubsection{概述}
三体问题是指空间中三个质点仅在相互的万有引力作用下运动的规律问题，现代数学告诉我们三体问题没有解析解，仅在限定条件的情况下（比如将三个质点限制在同一平面内，并且使其中一个质点的质量远小于另外两个质点），通过近似算法能够得到三体问题的一些特解。自从1900年数学家希尔伯特提出N体问题的特例——三体问题，一百多年来，不断有数学家和物理学家找出三体问题的特殊解，而在所有的解中，最重要也最常见的莫过于平面圆形限制性三体问题的五个特解，即拉格朗日点，无论是在古老的天体力学理论中，还是在如今前沿的太空探测研究中，拉格朗日点都具有弥足珍贵的重要地位。

本文将通过分析拉格朗日系统的势能分布情况，先找到5个拉格朗日点的大致位置，再给出拉格朗日点位置表达式的近似推导方法。
\subsubsection{两体近似}
上文提到三体问题没有解析解，即便是在平面圆形限制性三体问题的约束条件下，给出拉格朗日点位置的解析表达式仍然是不可能的，由此产生了两体近似的方法。两体近似是完全忽略质点3对质点1与质点2的影响（3个质点质量满足$m_1\gg m_2\gg m_3$），质点1、2在万有引力作用下绕两者质心O做匀速圆周运动，距离为R，角速度为$\omega=\dot{\theta}$，引入二体约化质量后的运动方程如下

\begin{align}
	\frac{Gm_1m_2}{R^2}&=\frac{m_1m_2}{m_1+m_2}R\dot{\theta}^2\label{LagrangePointsEQ001}\\
	\frac{G(m_1+m_2)}{R^3}&=\dot{\theta}^2\label{LagrangePointsEQ002}
\end{align}

解得二体转动角速度为。由于质点3质量相对于质点1、2而言非常小，我们有理由近似地认为，当质点3被引入后，质点1、2的运动状态不变，即、R以及质心位置都不变。以下的分析及推导统一以二体近似为前提。
\subsubsection{拉格朗日系统的势能分布情况}
在相对质点1、2静止的转动参考系下，拉格朗日体系的势场由两部分构成，其一是质点1、2自身产生的引力势，其二是转动系下的惯性离心势。为了定量地考察两种势能，以质心为坐标原点，以质点1、2连线为x轴建立平面直角坐标系，质点2在x轴正半轴上。

\begin{figure}
	\includegraphics[scale=0.3]{2BodiesSystem}
	\caption{建立坐标系\label{2BodiesSystem}}
\end{figure}

如图\ref{2BodiesSystem}所示，以O为势能零点，记惯性离心势为$U_{cf}(x,y)$，引力势为$U_G(x,y)$，设任意一点p(x,y)与质点1、质点2的距离分别为$l_1,l_2$，由几何关系得

\begin{align}
	l_1&=\sqrt{y^2+(\frac{m_2}{m_1+m_2}R+x)^2}\label{LagrangePointsEQ201}\\
	l_2&=\sqrt{y^2+(\frac{m_1}{m_1+m_2}R+x)^2}\label{LagrangePointsEQ202}
\end{align}

而惯性离心势与引力势分别为

\begin{align}
	U_{cf}(x,y)&=\int_{\sqrt{x^2+y^2}}^{0}\dot{\theta}^2rdr\\
	&=-\frac{G(m_1+m_2)}{2R^3}(x^2+y^2)\label{LagrangePointsEQ203}\\
	U_G(x,y)&=-\frac{Gm_1}{l_1}-\frac{Gm_2}{l_2}+\frac{Gm_1}{\frac{m_2}{m_1+m_2}R}+\frac{Gm_2}{\frac{m_1}{m_1+m_2}R}\label{LagrangePointsEQ204}
\end{align}

势场分布为

\begin{align}
	U(x,y)&=-\frac{G(m_1+m_2)}{2R^3}(x^2+y^2)-\frac{Gm_1}{l_1}-\frac{Gm_2}{l_2}+\frac{Gm_1}{\frac{m_2}{m_1+m_2}R}+\frac{Gm_2}{\frac{m_1}{m_1+m_2}R}\label{LagrangePointsEQ205}
\end{align}

为方便对势函数进行分析，定义无量纲变量如下：

\begin{align}
	\Phi&=\frac{RU}{Gm_1}\label{LagrangePointsEQ206}\\
	X&=\frac{x}{R}\label{LagrangePointsEQ207}\\
	Y&=\frac{y}{R}\label{LagrangePointsEQ208}\\
	\alpha&=\frac{m_2}{m_1}\label{LagrangePointsEQ209}
\end{align}

以上四式带入\ref{LagrangePointsEQ205}式并化简得

\begin{align}
	\Phi(x,y)&=-\frac{1+\alpha}{2}(X^2+Y^2)-\frac{1}{\sqrt{Y^2+(\frac{\alpha}{1+\alpha}+X)^2}}-\frac{\alpha}{\sqrt{Y^2+(\frac{1}{1+\alpha}-X)^2}}+\frac{1}{\alpha}+1+\alpha+\alpha^2\label{LagrangePointsEQ210}
\end{align}

利用MATLAB做出三维等势线图\ref{3DisopotentialMap}与二维等势线图\ref{2DisopotentialMap}

\begin{figure}
	\includegraphics[scale=0.3]{3DisopotentialMap}
	\caption{三维等势线图\label{3DisopotentialMap}}
\end{figure}

\begin{figure}
	\includegraphics[scale=0.3]{2DisopotentialMap}
	\caption{二维等势线图\label{2DisopotentialMap}}
\end{figure}

观察等势线图发现，在质点1、2之间，质点2的右边，质点1的左边以及质点1、2连线的上下两侧，势分布呈现出“鞍”、“峰”或“谷”的特性，也就是说在这些位置附近存在质点3的平衡位置，也即拉格朗日点。

\begin{figure}
	\includegraphics[scale=0.3]{LagrangePointsPos}
	\caption{拉格朗日点的位置\label{LagrangePointsPos}}
\end{figure}

在质点2的左右两侧，并靠近质点2的位置分别是拉格朗日点L1和L2，设其相对于质点2的距离分别为r1和r2（待求），在质点1的左侧靠近点2'的位置是L3（点2'与点2关于点1对称，这样设定是为了方便后续计算与表达），设其相对于点2'的距离为r3（待求）。在质点1、2的上方（第一象限）与下方（第四象限）分别为拉格朗日点L4和L5。
\subsubsection{计算L1、L2、L3的位置}
令

\begin{align}
	\beta_1&=\frac{r_1}{R}\label{LagrangePointsEQ214}
\end{align}

质点3在L1处的力学方程如下

\begin{align}
	m_3\dot{\theta}^2(\frac{m_1}{m_1+m_2}R-r_1)&=\frac{Gm_1m_3}{(R-r_1)^2}-\frac{Gm_2m_3}{r_1^2}\label{LagrangePointsEQ211}\\
	\dot{\theta}^2R(\frac{1}{1+m_2/m_1}-r_1/R)&=\frac{Gm_1}{R^2}(\frac{1}{(1-r_1/R)^2}-\frac{m_2/m_1}{(r_1/R)^2})\label{LagrangePointsEQ212}
\end{align}

式\ref{LagrangePointsEQ209}、\ref{LagrangePointsEQ214}、\ref{LagrangePointsEQ002}代入\ref{LagrangePointsEQ212}，得到

\begin{align}
	\frac{G(m_1+m_2)}{R^3}(\frac{1}{1+\alpha}-\beta_1)&=\frac{Gm_1}{R^3}(\frac{1}{(1-\beta_1)^2}-\frac{\alpha}{\beta_1^2})\label{LagrangePointsEQ216}\\
	(1+m_2/m1)(\frac{1}{1+\alpha}-\beta_1)&=\frac{1}{(1-\beta_1)^2}-\frac{\alpha}{\beta_1^2}\label{LagrangePointsEQ217}\\	
	1-(1+\alpha)\beta_1&=\frac{1}{(1-\beta_1)^2}-\frac{\alpha}{\beta_1^2}\label{LagrangePointsEQ218}\\	
	\beta_1^2-(1+\alpha)\beta_1^3&=\frac{\beta_1^2}{(1-\beta_1)^2}-\alpha\label{LagrangePointsEQ219}
\end{align}

注意到$\alpha\ll1,\beta_1\ll1$，可以对上式等号右边的第一项进行近似计算，之后约去高阶小量得到一个近似结果

\begin{align}
	\beta_1^2-(1+\alpha)\beta_1^3&=\beta_1^2(1+2\beta_1)-\alpha\label{LagrangePointsEQ220}\\
	\alpha&=\beta_1^3(3+\alpha)\label{LagrangePointsEQ221}\\
	\alpha&=3\beta_1^3\label{LagrangePointsEQ222}\\
	\beta_1&=\sqrt[3]{\frac{\alpha}{3}}\label{LagrangePointsEQ223}	
\end{align}

$\alpha=3\beta_1^3$，即$r_1=R(m_2/m_1)^{1/3}$（这一结果满足了$r1\ll R$，说明我们的近似没有过度，同时还发现$\alpha\sim\beta_1^3$这一特点）。

L2与L1完全类似，这里不再过多赘述，直接给出结果，，即，结果说明与的位置与质点2对称！
质点3在处的力学方程与前两个方程的差异性较大：

仍然通过近似计算解的，即（这一结果与前两个差异较大，但仍满足，但是）。
、的位置推导
在推导和的位置之前，先考虑更一般的情况：平面内任意三个质点成三角形绕其公共质心旋转并保持相对静止，如图所示，我们先讨论这样的稳定体系的几何结构。

三角形体系
质点1受到质点2、3的引力之和，即

质心位置公式

又因为质点1绕质心c做匀速圆周运动，故质点1所受的引力之和应当指向质心，即，故有如下关系式

化简得，同理可得，故可知质点1、2、3之间构成等边三角形！
显然和即为上述结的特殊情况（），故和各自与质点1、2构成等边三角形。
\subsubsection{小结}
本文首先通过分析拉格朗日点的势场分布情况定性地确定了五个拉格朗日点的位置，再通过力学方程的近似求解解得五个拉格朗日点的具体位置，在上述坐标系下，拉格朗日点的坐标表示如下：





发布于2019-04-06
